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ABSTRACT 

Meridional flows with velocities of a few meters per second are observed in the uppermost regions of the solar 
convection zone. The amplitude and pattern of the flows deeper in the solar interior, in particular near the top 
of the radiative region, are of crucial importance to a wide range of solar magnetohydrodynamical processes. In 
this paper, we provide a systematic study of the penetration of large-scale meridional flows from the convection 
zone into the radiative zone. In particular, we study the effects of the assumed boundary conditions applied at 
the convective-radiative interface on the deeper flows. Using simplified analytical models in conjunction with 
more complete numerical methods, we show that penetration of the convectively-driven meridional flows into the 
deeper interior is not necessarily limited to a shallow Ekman depth but can penetrate much deeper, depending on 
how the convective-radiative interface flows are modeled. 

Subject headings: hydrodynamics — method:numerical — method:analytical — Sun:interior 



1. INTRODUCTION 

Meridional flows in the solar interior have recently become 
the focus of observational and theoretical attention. Poleward 
sub-surface flows with amplitudes of the order of a few tens 
of meters per second have been detected with reliable accuracy 
down to about O.85r (Giles et al. 1997). A globally equator- 
ward return flow must exist deeper in the interior to guarantee 
mass conservation, but its amplitude and structure can only be 
conjectured currently. This paper addresses the question of how 
deeply these meridional flows penetrate into the radiative zone. 
An understanding of these return meridional flows is of fun- 
damental importance since their nature plays a crucial role in 
many current theories for the internal magneto-hydrodynamics 
of the Sun. 

Firstly, meridional circulations have been argued to play a 
central role in the operation of the global solar dynamo (see the 
review by Charbonneau, 2005). In these models, the predicted 
spatio-temporal behavior of the solar cycle depends sensitively 
on the assumed circulation pattern and speed. The chosen posi- 
tion for the return flow coupled with mass conservation sets the 
velocity of the equatorward flow near the base of the convec- 
tion zone and therefore controls the activity cycle period. Sim- 
ilarly, the depth of penetration of the meridional flows into the 
radiative region determines where the toroidal magnetic field 
is generated, and therefore also influences the cycle period and 
the field amplitudes. 

Secondly, meridional flows advect angular momentum, and 
therefore play a key role in the global dynamical balance of the 
solar interior. For example, helioseismology has revealed the 
existence of a strong radial shear layer, now known as the solar 
tachocline (Brown et al. 1989; Spiegel & Zahn, 1992; Hughes, 
Rosner & Weiss, 2007), located precisely at the interface be- 
tween the radiative and convective regions. Quantitative mod- 
els of the tachocline have revealed a sensitive dependence of 
the interior angular velocity profile on the derived or assumed 
interfacial flows (Spiegel & Zahn, 1992; Gough & Mclntyre, 
1998; Rempel, 2005; Garaud, 2007). 

Finally, meridional flows also transport various chemical 
species within the solar interior, with directly observable con- 



sequences. Near the base of the convection zone, mixing by 
large-scale flows can prevent the gravitational settling of helium 
with respect to hydrogen, leaving a noticeable signature in the 
helioseismic sound-speed data (Elliott & Gough ,1999). Addi- 
tional mixing below the convection zone in the main-sequence 
phase is also required by the observed surface abundances of 
light elements such as lithium and beryllium, with plausibly the 
same origin. 

Flows with amplitudes on the order of tens of meters per sec- 
ond are required to balance angular momentum transport by the 
turbulent stresses throughout the convection zone (Miesch et al. 
2000). Since there exist no physical barrier between the con- 
vective and radiative regions, these convectively driven flows 
may continue their downward progress somewhat beyond the 
driving region into the radiative zone, thereby "penetrating" or 
"overshooting" into the interior while retaining potentially sig- 
nificant velocities. How far flows extend beyond their driving 
region is an essential question. 

The problem has recently been addressed by Gilman & Mi- 
esch (2004) (GM04 hereafter) and by Mclntyre (2007). Using a 
steady-state formalism, GM04 discuss the depth of penetration 
of an existing latitudinal flow into the radiative interior. They 
argue that the source flow amplitude is rapidly damped in the 
radiative interior within a shallow Ekman depth. This depth 
ranges from a fraction of a kilometer, using microscopic values 
of the viscosity, to a few tens of kilometers, using a turbulent 
value of the viscosity. Gilman and Miesch reach a strong con- 
clusion, namely that "the physics of the solar tachocline and 
neighboring regions does not allow penetration of meridional 
circulation originating in the solar convection zone below the 
overshoot layer". This could have dramatic consequences for 
the magneto-hydrodynamics of the solar radiative zone. 

However, steady-state solutions are sensitively dependent on 
boundary conditions. GM04 solve the problem for the radia- 
tive interior dynamics where the forcing is by purely latitudinal 
flows at the upper convective-radiative interface. They do not 
consider interfacial forcing from flows generated by the differ- 
ential rotation nor by direct radial pumping into the radiative 
zone (only radial flows generated for mass conservation in re- 
sponse to their latitudinal forcing). It is reasonable to address 
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the question as to whether their strong conclusion remains ap- 
plicable under more general circumstances. 

In this paper, we therefore extend the work of GM04 to al- 
low for greater generality in the source forcing flow, allowing 
the possibility of azimuthal and radial flows in addition to the 
latitudinal flows. We systematically examine the consequences 
of using these various sets of boundary conditions to mimic 
the convective-radiative interface. We find that the meridional 
flows can penetrate to significantly different depths, depending 
upon the choice of boundary conditions. The GM04 solution 
can be recovered in special cases, but is typically overpowered 
by other solutions when alternative assumptions are made about 
the nature of the flows driven within the convection zone. 

In what follows, we adopt a three-step approach to study the 
penetration of meridional flows into the solar radiative zone. 
Given the added difficulties inherent to spherical coordinate 
systems, we first examine the problem in Cartesian geometry. 
In §|2] we study analytically the complete set of steady, linear, 
pseudo-axisymmetric Boussinesq equations and explore a wide 
range of boundary conditions that may possibly mimic the ef- 
fects of the convection zone on the radiative zone. We system- 
atically discuss the solutions obtained, which are linear combi- 
nations of two fundamental modes of behavior. One of these 
modes is a solution which varies rapidly on a typical shallow 
Ekman scale as found by GM04, the other one is a more slowly 
varying solution which can span the entire interior. In order to 
gain better insight into the physics of the system and in partic- 
ular the new solutions, in §[3] we consider a simplified Carte- 
sian Boussinesq model in which the solutions related to the Ek- 
man layers are artificially suppressed by neglecting the viscous 
terms in the radial and latitudinal components of the momen- 
tum equation. Finally, in §0] we relax the Cartesian constraint 
and present numerical results for all of the various boundary 
conditions in a steady, linear, anelastic, spherical but axisym- 
metric simulation of the solar radiative zone. We compare these 
numerical results in spherical geometry with the analytical pre- 
dictions from the Cartesian models, both in terms of the scale of 
variation of the solutions, and in terms of their predicted flow 
velocities. While the spherical geometry as well as the non- 
uniform background state necessarily add to the complexity of 
the problem, we find that the analytical scalings extracted from 
the Cartesian geometry models agree very well with the full 
numerical solutions. When more complex boundary conditions 
are taken into account, limits on the penetration of meridional 
flows into the radiative zone are much less stringent than previ- 
ously claimed. The exact flow velocities and depth of penetra- 
tion achieved, however, depend sensitively on the actual bound- 
ary conditions selected. These results and conclusions are dis- 
cussed in detail in SB] 



2. A CARTESIAN MODEL 

In all that follows we consider a stably stratified radiative 
zone located beneath a turbulent convection zone, a typical sit- 
uation encountered in all solar- type stars. Within the convection 
zone, we assume that turbulent stresses drive large-scale flows 
in the azimuthal direction (i.e. a large-scale differential rota- 
tion) as well as in the meridional direction. The amplitude and 
spatial variation of the flows just above the convective-radiative 
interface is assumed to be known. We then pose and solve the 
following question: What is the resulting flow pattern and ve- 
locities in the underlying radiative zone? 



2.1, Model setup 

The spherical geometry of the solar radiative zone as well as 
the non-uniform background state (in terms of temperature and 
density, viscosity and thermal conductivity for instance) both 
preclude any attempt at solving the problem analytically. We 
postpone to §|4]the presentation and discussion of the complete 
numerical solution of the problem, and first consider a much 
simplified "radiative zone" with rectangular geometry in Carte- 
sian coordinates (x,y,z), and a uniform background temperature 
gradient, viscosity and thermal conductivity. In this coordinate 
system, the x-direction can be thought of as the azimuthal di- 
rection with x E [Q,2ttR], the y-direction is aligned with the 
latitudinal direction and is limited to y £ [0, irR] and finally the 
z-direction is the radial direction with z<R- The poles are rep- 
resented by y = and y = ttR while the equator is at y = wR/2, 
The dimensional constant R represents the base of the convec- 
tion zone, and z = the interior of the Sun. The system rotates 
with angular velocity O, = (0, 0, 0), and gravity is assumed to be 
aligned with the rotation axis. 

2.2. Model equations and general solution 

The equations governing dynamical and thermal perturba- 
tions to a stably stratified background assumed at rest are 
the mass, momentum and thermal energy conservation equa- 
tions. Using the Boussinesq approximation and assuming "ax- 
ial" symmetry (i.e. d/dx = 0), we first linearize these equations 
in the thermal perturbations and flow velocities, then project 
them onto the Cartesian coordinate system as 
dv + dw 

dy dz 

2il { ^ M + 

\ dy 2 dz 2 

20, ^ + ( + 

dy \ dy 2 dz 2 

= -^ + ga9, 
dz 

where u = (u,v,w) is the flow velocity, is the temperature 
perturbation, v and k are the viscosity and thermal diffusivity, 
a = l/T where T is the background temperature and finally, 
ga(3 = N 2 is the background buoyancy frequency. Note that the 
viscous diffusion term in the radial component of the momen- 
tum equation has been neglected in accordance with hydrostatic 
equilibrium; this does not affect in any way the conclusions of 
this paper. 

While GM04 and Mclntyre (2007) neglected the d 2 /dy 2 
terms in the viscous terms, we consider here for completeness 
the full Laplacian. The additional terms are found to be nec- 
essary in the light of the fact that some of the boundary layers 
in the system are large compared with the vertical size of the 
domain. Neglecting these terms is not physically justified, and 
mathematically transforms any slowly varying, standard expo- 
nential solutions into the rather unphysical secular linear solu- 
tions described by Mclntyre (2007) (and neglected by GM04). 

The hydrostatic equilibrium equation can be combined with 
the latitudinal component of the momentum equation to yield a 
generalized thermal-wind equation, 

2 n — + a — -v— (— + — \ (2) 

dz aS dy dz \dy 2 dz 2 J 
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Seeking solutions with latitudinal dependence as sin(2ny/R) 
or cos(2ny/R), and exponential vertical dependence as exp(fe), 
we obtain the characteristic equation 
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where we have introduced two standard characteristic length- 
scales 
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andJ BD =(p) 
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the first one being the standard Ekman depth and the second 
representing a buoyancy-diffusion layer (GM04, Barcildon & 
Pedlosky 1967). 

It is possible to show (with some algebra) that in the limit 
where 
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(which is always true below the base of the convection zone) 
then the eight solutions to equation (01 (±&i, ±&2> ±£3 and ±£4) 
can be approximated by 

k\ = — , 
R 

dl 2n 
h ~ -f , 

4d r 

h A ~^{\±i)d£, (6) 

thus yielding four slowly varying exponential solutions (related 
to k\ and k 2 ) in addition to four rapidly varying, oscillatory and 
exponential solutions related to £3 and k 4 which describe the 
Ekman layers of the system. The Ekman solutions were found 
and described by GM04. On the other hand, the ±k\ and ±£2 
solutions differ from the (two) real solutions found by GM04, 
a discrepancy which can easily be traced back to the omitted 
d 2 /dy 2 in their viscous diffusion terms. 
It is interesting to note that 

k 2 = VPrBu^ = VPrBu^k! , (7) 

where D is the local density scaleheight, Pr = v/k is the Prandtl 
number and where the Burger number Bu is defined as 
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The Burger number in the solar radiative zone is estimated to be 



about Bu : 



2.5 x 10 3 using = 2.7x 10^s" 



,JV = 8 x KTV 



R = 5 x 10 10 cm and D = OMR = 8.6 x 10 9 cm (see Gough, 
2007). In laminar regions of the solar radiative zone, the mi- 
croscopic Prandtl number is of order Pr = 2x 10~ 6 , so that 

k~ 2 l ~ 5k\ x , (9) 
which is clearly of the order of R itself for large-scale forcing 
(small n). 

2.3. General solutions 

Consider for instance the solution of (03 corresponding to 
cos(2ny/R) and sin(2ny/R) variations in the azimuthal veloc- 
ity: 



u(x,y,z)= (a\e k,z + a 2 e kiZ + a 3 e klZ + a 4 e klZ 



+ a 5 e hz +a 6 e hz + a 1 e hz + a%e hz ) cos 



(b ie klZ + b 2 e~ klZ + b 3 e k2Z + b 4 e~ klZ 



(10) 
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This corresponds to 

v(x,y,z) = dl [(ki-kf) ( ai e k2Z + a 4 e- k2Z ) 
+ (k 2 -kl) (a 5 e hz + a 6 e- k *) 

+ (k 2 -kf) (a 7 e k4Z + a & e~ k4Z )] cos 

+ d\ [(k\-k 2 2 ) (b 3 e k2Z + b 4 e- k2Z ) 
+ (k 2 -k 2 ) (b 5 e k ' z + b 6 e- k ' z ) 

+ (k 2 -kj) (b 7 e^+he-^)] S m(^\ 
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w(x,y,z) = d E 
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K3 
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K 4 
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and finally 

0(x,y,z) = 
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where one must bear in mind that k[ and k 2 depend on the 
wavenumber n/R of the forcing function. 

2.4. Boundary conditions 

To find the amplitude of each exponential term, one must ap- 
ply the boundary conditions to the general solutions. A quick 
look at the system shows that eight boundary conditions are re- 
quired, which arise from conditions on u, v, w and 9 at both the 
top and bottom of the domain. 

We choose boundary conditions at z = R to represent the ac- 
tion of the convection zone on the underlying stably stratified 
and laminar radiative region. We require the continuity of the 
radiative interior solution with the complete vector of velocities 
at the base of the convection zone, so that 



u(x,y,R) = u cz (y) = (u cz (y), v cz (y),w cz (y)) 



(14) 
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(where the meridional and azimuthal flows in the convection 
zone are assumed to be axisymmetric). A simple reasonable 
prescription for the flow velocities at the interface might be, for 
instance, 



u cz (y) = -U cos 



Vcz(y) = V sin 



2v 



w cz (y) = -W cos 



(15) 



The azimuthal forcing term u cz (y) represents a solar-like differ- 
ential rotation (with slower-than-average rotation near the poles 
and faster-than-average rotation near the equator if Uq > 0). The 
latitudinal and radial forcing terms v cz (y) and w cz (y) represent 
a single-cell flow in each hemisphere with independent verti- 
cal and latitudinal flow velocities. When Vb > the flow is 
equatorward at the boundary, while Wo > guarantees inflow 
in the high latitudes and outflow in the low latitudes. Note that 
since the problem studied is linear, and since any of the three 
velocity components of (u cz (y), v cz (y),w cz (y)) can be written as 
a Fourier series in cos(2ny/R) and sin(2ny/R), it is possible to 
find the general solution of (Q~|i for any set of imposed velocities. 
The profiles chosen here contain only one Fourier component 
for simplicity. 

Near the bottom boundary, we would in principle like to 
choose boundary conditions that have as little effect on the solu- 
tion as possible. In the real Sun they would instead be replaced 
by regularity conditions at the origin. However, fitting eight 
boundary conditions to the general solutions yields an 8 x 8 lin- 
ear system which is difficult to study analytically (even with the 
help of Maple or equivalent software). Thus in this section we 
restrict our study to the case where the bottom boundary is lo- 
cated at z — > — oo. This immediately implies that all of the even 
{a,} and {/?,} coefficients must be null. In $3] we revisit this 
simplification, and suggest an easy way of deducing the solu- 
tions in a geometry where the bottom boundary is at z = from 
those in a semi-infinite domain. In $4]we verify our hypothesis 
against numerical results. 

In what follows ( jj2.5| and 32. 6t we describe two possibilities 
for the thermal boundary conditions at the convective-radiative 
interface. 

2.5. Thermal boundary conditions of type 1: "perfectly" 
conducting convection zone 

The local heat flux through the boundary associated with the 
perturbations is the sum of the conducted heat flux -kd6/dz 
and the advected heat flux phw, where p is the background den- 
sity, h is the background enthalpy, and k = ~pc p K is the thermal 
conductivity (here c p is the specific heat at constant pressure). 

In a steady-state this local heat flux must be equal to zero, so 
that when the convection zone is perfectly conducting (k — > oo) 
the thermal boundary condition reduces to dO/dz = 0. 

Using this final boundary condition together with the ones 
described in 32.41 we find that the remaining odd coefficients 
{<*i}i=i, 3.5,7 satisfy the linear system MA = C with 
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and 
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Similarly, the coefficients {£>/},=i.3,5.7 satisfy the systemMZ? = D 
with 

fbie klR \ / 
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In the limit where ki -C k\ -C |&3|, \k^\ the solutions can be 
simplified to yield 
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and 
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We can now finally evaluate, for instance, the latitudinal flow 
velocity within this rectangular radiative zone: 



v(x,y,z) 



-dlk\Uoe kl(z 



-R). 



2-2^0 k A k 3 Kfe-/Q 
El' ■ ■ 



d\k\ 



Up hk A 
£2 £4 — ^3 
k 2 

-Vot^- + Wo 

K3K4 



Mz-R) 



V, 



w 



ki ki-k4 
v R 

k\ (k-i + k$) 
k^4 



(21) 
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As expected from the linearity of the governing system equa- 
tions (fl]i, the meridional flow velocity in the radiative interior 
is simply the sum of the three contributions arising from az- 
imuthal forcing only (with terms proportional to Uo), latitudinal 
forcing only (with terms proportional to Vo) and radial forcing 
only (with terms proportional to Wo). In addition, each of these 
three contributions is the sum of three terms, one with exponen- 
tial dependence in e kl( - z ~ R) which corresponds to a very slowly 
varying function of depth, and two complex conjugate terms 
with exponential dependence in e kl( - z ~ R ^ and e ki( - z ~ R) associated 
with the very rapidly decaying and oscillating Ekman solutions. 

We now compare the relative amplitudes of all of these terms 
as a function of the parameters of the system. The amplitude 
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of the rapidly decaying component of the solution arising from 
azimuthal forcing only is 



d 2 k 2 ^ 
E 1 k 2 



k 3 k 4 



k 3 — k 4 



8ff„ D 
PrBu R 



:U 



(22) 



using the values of ki derived in equations © and (0, and 
where 



_ . _2dl 



v 



(23) 



is the Ekman number. Similarly, the slowly decaying compo- 
nent of the solution arising from azimuthal forcing only has an 
amplitude 

4kjU ~ 2E V U (24) 

Given that E v ~ 10~ 16 near the base of the convection zone for 
microscopic values of the viscosity, within an Ekman length of 
the boundary both components of the solution are negligible. 

Similarly with the other forcing contributions, we conclude 
that 

• latitudinal forcing drives meridional flows with a rapidly 
decaying component which has an amplitude Vb (with 
no dependence on any of the other parameters of the sys- 
tem), while the slowly decaying component has an am- 
plitude proportional to E u Vq. Thus, in agreement with 
the study of GM04, we find that within a few Ekman 
lengths, both are negligible for microscopic values of 
the viscosity. 

• radial forcing drives meridional flows with a rapidly de- 
caying component which has an amplitude proportional 
to Wo/t/E^, while the slowly decaying component has 
an amplitude proportional to ^/E^Wq. Thus, in this par- 
ticular case we find that beyond a few Ekman lengths 
the slow mode retains a non-negligible amplitude of 
about 10~ 8 times the imposed velocity. Indeed, for im- 
posed meter per second flow velocities, the flows near 
the top of the radiative zone could then be of the or- 
der of i/EZWq ~ 10~ 6 centimeters per second, with an 
overall turnover time of the order of R/ t/E^Wq ~ 1 Gyr. 
While very slow, this can still provide mixing in the 
tachocline on the stellar evolution and/or gravitational 
settling timescale. 

2.6. Thermal boundary conditions of type 2: no net perturbed 
heat flux through the boundary 

Dropping the assumption that the convection zone is a per- 
fectly conducting fluid, we require instead that hw = c v ndO/dz. 
which is equivalent to Tw = ndd/dz. 

We can rewrite this new thermal boundary condition into the 
linear systems MA = C and MB = E where M, A, B and C have 
already been defined, and 
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where H® is the potential temperature scaleheight 
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Thus the {a,} coefficients are the same as in the previous sec- 
tion, while in the limit where kj <C k\ <C \ks\, \k 4 \ 
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For this new set of boundary conditions, azimuthal forcing 
and latitudinal forcing yield the same solutions as in the previ- 
ous section, while radial forcing drives meridional flows with 
a rapidly decaying component which has an amplitude propor- 
tional to Wo(Hsk 2 - 1)/ \[eZ, while the slowly decaying compo- 
nent has an amplitude proportional to {Hq / R){R / D) 2 PrBuWo- 
Note how in this case the slowly decaying component of the 
flow has an amplitude which is independent of the background 
viscosity v. 

2.7. Discussion of the solutions 

Using this Cartesian geometry model, we have shown that 
there exist solutions for meridional flows in the radiative in- 
terior with significant amplitude throughout. These flows can 
only be present if driven by direct pumping into and out of the 
radiative zone, that is, when w cz (y) has a significant amplitude 
on the boundary. Our study also shows that the actual ampli- 
tude of the flows within the radiative zone depends sensitively 
on the thermal boundary conditions used. While it is not clear 
which set of boundary conditions actually accurately represents 
the true convective-radiative interface, it is not implausible that 
they may lie somewhere in between the two extreme cases stud- 
ied. Thus it is also not implausible that there may be significant 
penetration of the convective zone flows into the radiative zone, 
with typical turnover times shorter than R/^E^Wo- 

Our results extend the work of GM04, and naturally recover 
their solutions in the limit where only latitudinal forcing is 
taken into account. Note that GM04 do not specifically re- 
quire that w be null on the convective-radiative interface, and 
argue that any spatially varying latitudinal flow drives radial 
flows into and out of the boundary to guarantee mass conser- 
vation. This is indeed correct, but one must keep in mind that 
dv/dy = -dw/dz only relates latitudinal variations in v to the 
radial derivative of w. Since the derivative is dominated at the 
base of the convection zone by rapid variation on the Ekman 
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scale, the mass conservation equation in fact requires that the 
radial flows in the GM04 solution have an amplitude of the or- 
der of Wo = VqcIe/R. It is therefore easy to see why their solution 
yields negligible velocities for the slowly varying component of 
the flow penetrating into the radiative zone. 

Finally, it is crucial to note that since the momentum and en- 
ergy equations in this study (and the preceding ones by GM04 
and Mclntyre, 2007) have been linearized, the predicted am- 
plitudes of the flows are linearly dependent on the imposed 
flow velocities. In reality, this will only be true in the limit 
where the typical amplitudes of the nonlinear advection terms 
(u ■ Vu for the momentum equation, and u ■ VO for the energy 
equation) are indeed much smaller than the corresponding lin- 
ear terms (2VL x u for the momentum equation and (3w for the 
energy equation). The nonlinear terms in the momentum equa- 
tion for example guarantee that no unphysical counter-rotation 
is allowed. The linearized equations on the other hand scale 
arbitrarily with the imposed boundary conditions and do allow 
counter-rotation for certain input parameters, e.g. if the im- 
posed differential rotation profile u cz (y) is large enough or if 
the meridional flows are too rapid (see Figure [5] for instance). 
Similarly, the nonlinear terms in the energy equation guaran- 
tee that the actual turnover time of the meridional flows cannot 
be lower than the local thermal diffusion time (over the depth 
considered), while the linearized equations allow for any val- 
ues of the turnover time provided the imposed boundary flow 
velocities are high enough. 

From this linearized model, we can therefore say with rea- 
sonable confidence that meridional flows can indeed penetrate 
into the radiative zone, provided their turnover time (typically 
estimated as v(d) jd where d = R-z is the depth considered) 
is longer than the local thermal diffusion time (estimated as 
d 2 /n). 



3. A REDUCED CARTESIAN MODEL 

One of the remaining issues that needs to be addressed is 
that of lower boundary. In the previous section, the high- 
dimensionality of the solution space made it difficult to con- 
sider the more realistic situation of a bottom boundary located 
at z = 0. While the Ekman solutions (associated with k 3 and 
k 4 ) decay inward/downward so quickly that the presence or ab- 
sence of a lower boundary cannot affect them, the slowly vary- 
ing solutions do span the entire radiative region from z = R down 
to z = 0; the location of the lower boundary is expected to have 
some influence on their amplitude. 

Here, we therefore focus on studying the slowly varying so- 
lutions only by considering a system in which the Ekman flows 
are filtered out. This is another way of reducing the dimen- 
sionality of the solution space, and enables us to compare the 
solutions in a semi-infinite domain to the solutions in a finite 
domain. 



We consider the system 



dv ^ dw 
dy dz 

dp 

& 



2Qu ■ 



dp 

dy 



-2Q.V-- 



fiw = K 



d 2 u d 2 u 

dy 2 dz 2 

d 2 9 d 2 6 

dy 2 dz 2 



(29) 



in which the viscous diffusion terms within the latitudinal com- 
ponent of the momentum equation has now been removed. This 
simplification is consistent with the assumptions of geostrophic 
equilibrium and, as we now prove, effectively filters out the 
Ekman flows. The system can be reduced to a single partial 
differential equation for u, for instance, as 



d 2 u di d A u 



dz 2 dy 2 d A m dy 4 



■ 0. 



(30) 



Seeking exponential solutions in z with periodic behavior in 
y as before yields the characteristic equation: 

4 „2\ / A „2 



R 2 



An 1 d\ 



■ . 



(31) 



with solutions ±K\ and ±K 2 with 

In 



K\ = — = k] and 
R 



K 2 = 



In 



R \d\ 



d l 



= k 



2 ■ 



(32) 



Thus we recover only the slowly varying solutions of the previ- 
ous section. 

3.2. General solutions 

The flow solution to the above system for fixed latitudinal 
wavenumber n is 



u(x,y,z) = (Aie KiZ +A2e- KiZ +A 3 e K2Z +A 4 e~ K2Z )cos 



2ny 
~R~ 



+ (B 1 e KiZ + B2e~ KiZ + B } e Klz + B 4 e- K2Z )sm ^^^33) 
which also yields 

v(x,y,z) = d 2 E (K 2 -K 2 ) (A,e K2Z +A 4 e- K2Z ) cos f ^ 



+ dl (K 2 -K 2 ) (B 3 e K2Z + B 4 e- K2Z ) sin (34) 



and 



,K, 



w(x,y,z) = dl 1 ^- <^ 2 -^ (a^ zZ _ a _„-k^ 



(Kt-Ki) (A 3 e ^-'-A 4 e~ 



sin 



2ny 
~R~ 



dlf 2 {Kl-K 2 ) [B.^ /> V * )cos (^35) 



3.1. Reduced model equations 



3.3. Boundary conditions 

The governing equations considered form what is apparently 
a 6 th order system in the z direction, and require, in principle, 
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6 boundary conditions split between the top and bottom bound- 
aries. Naturally, one would like to impose boundary conditions 
on u, w and 9 (one at the top, one at the bottom for each vari- 
able). Note that since the viscous terms have been neglected 
in the equations for hydrostatic and geostrophic equilibrium re- 
spectively, one cannot impose a boundary condition on v: the 
equations contain no stresses that could transfer the boundary 
information to the rest of the system. However, we see that 
combining these equations only yields a 4 th order partial differ- 
ential equation for u. This implies that the 6 selected boundary 
conditions must somehow be redundant, otherwise there will be 
no solution to the system. In what follows, we therefore only 
select boundary conditions for u and w at the top and bottom 
boundaries. Near the top boundary, we consider as before 

u(x,y,R) = u cz (y) , 

w(x,y,R) = w cz (y) . (36) 

When considering a bottom boundary at z = 0, we choose im- 
permeable boundary conditions for w, and stress-free boundary 
conditions for u so that 

du 



implying 



dz 



0. 



(x,y,0) 
w(x,y,0) = 0. 

3.4. Solutions for a semi-infinite domain 
When the bottom boundary is near — oo, we find that 
A, =-U e- KlR , 
A 2 =A^=A 4 = 0. 
Ki Wo KiR 



(37) 



B\ =- 



Br 



K\(Kf 



K 2 



-Kl) 4 

B 2 = B 4 
W _ 



■ 



K 2 R 



K\(Kf 



-Kl) 4 



so that 



K 2 

v(x,y,z) = —Woe 



K 2 (z-R) 



sin^y) , 



(38) 



(39) 



which illustrates again how radial forcing can yield non-zero 
flow amplitudes penetrating deeply into the radiative zone. It is 
important to note, however, that the predicted flow amplitude is 
different from that found in §|2] This might be attributed to the 
fact that when viscous effects are taken into account, a fraction 
of the flow penetrating into the radiative zone is deflected into 
the very shallow Ekman layers. 

3.5. Solutions for a finite domain 
When the bottom boundary is located at z = 0, we find that 

Uo 



At =A 2 



2cosh(KiR) 



A 3 = A 4 = 0, 



Bi = B 2 ■ 



K 2 Wo cosh(K 2 R) 
Ki (Kl - K\) ~4 2 smh(K 2 R) cos\\{K x R) 



Wo 



2dlK} KxRcosUKxR) 



B-k = Ba = 



K 2 Wo 



Kx{K\-Kl) d\ 2sinhCK 2 fl) 



Wq 1 
l4K\ KiR 



(40) 



v(x,y,z) = 



Wq K 2 
sinh(K 2 R) #i 



cosh(K 2 z) sin(Kiy) 



Wq 

— — cosh(K 2 z) sm( Kxy) . 



(41) 



3.6. Consequences 

Comparing the expressions in equations ( |39i > and ( PiTl i, we 
find that the slowly varying component of the meridional flow 
in the case of a finite domain has an amplitude that is 1 /K 2 R 
times that of the semi-infinite domain case. The difference be- 
tween the predicted amplitudes for the two geometrical systems 
(finite and semi-infinite domains) can easily be understood in 
the light of the fact that the exponential solutions associated 
with Ki and K 2 decay on much longer lengthscales than R. 

Extrapolating this result to the full Cartesian problem stud- 
ied in $2] we therefore predict that the slowly varying com- 
ponents of the meridional flows (associated with the k\ and k 2 
wavenumbers) should in fact have an amplitude that is 1 / k 2 R 
times that given in d22l i and ( |29l ; the Ekman components on 
the other hand decay so rapidly that their amplitude should not 
be influenced by the presence of a lower boundary. We now re- 
vise our estimates of 32.5| and j j2.6| to predict that in the case of 
type 1 boundary conditions (88 /dz = at the upper boundary) 
then 

• azimuthal forcing leads to meridional flows with a 
rapidly decaying component with an amplitude propor- 
tional to y / E u /PrBuUo, while the slowly decaying com- 
ponent has an amplitude proportional to E v Uo / \JPrBu. 

• latitudinal forcing leads to meridional flows with a 
rapidly decaying component with an amplitude propor- 
tional to Vb, while the slowly decaying component has 
an amplitude proportional to E v Vo / \fPrBu. 

• radial forcing leads to meridional flows with a rapidly 
decaying component with an amplitude proportional to 
Wo J \/E^, while the slowly decaying component has an 
amplitude proportional to y/E^/PrBuWo- 

In the case of type 2 boundary conditions (where ndd/dz = Tw 
at the upper boundary) the flow velocities predicted in the case 
of azimuthal and latitudinal forcing are the same, while radial 
forcing drives flows with a rapidly decaying component with 
an amplitude proportional to Wo(He>k 2 - \ )j\[EZ and a slowly 
decaying component which has an amplitude proportional to 
(H /R)VP?MWo. 



4. NUMERICAL SOLUTIONS IN A SPHERICAL SHELL 

4. 1 . Model setup 

To complete this axisymmetric study and illustrate our an- 
alytical results numerically, we perform simulations of flows 
in the solar radiative zone subject to various boundary condi- 
tions. The governing equations for the numerical model are 
derived by perturbing the spherically symmetric equations of 
stellar structure, moving to a rotating frame of reference and 
assuming that the velocity perturbations and thermodynamic al 
perturbations to the background spherically symmetric state are 
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small enough for linearization to be appropriate. Then 
2pTi x u = -V/J-pVl+ZV ■ n , 
V • (pu) = , 
pTu-Vi = /V-(Wf), 



where u = (u r ,ug,u^,) is the velocity field in a frame rotating 
with angular velocity ft, p, T, s and p are the standard thermo- 
dynamical variables, k = ~pc v n is the thermal conductivity, II is 
the viscous stress tensor (which depends on the viscosity v) and 
$ is the gravitational potential. Quantities denoted with bars 
are background quantities, taken from the standard solar model 
of Christensen-Dalsgaard et al. (1991), while the temperature, 
pressure and density perturbations are denotes with tildes. 

This system is fully consistent with the Cartesian model 
equations presented in with the added sophistication of the 
perfect gas equation of state and the anelastic approximation 
instead of the simpler Boussinesq approximation. This modifi- 
cation is added so that the model equations are consistent with 
the level of approximation used, but does not affect the nature of 
the solutions. The centrifugal force associated with the rotation 
of the background state fi x ft x r has been removed to suppress 
global Eddington-Sweet circulations, which are known to be of 
very small amplitude, but would otherwise play an important 
role in this steady-state calculation (see the work of Garaud, 
2002, for comparison). Perturbations in the gravitational poten- 
tial are neglected in accordance with Cowling's approximation. 

Note how both diffusion terms have been multiplied by the 
same factor /. Since the typical Ekman number in the Sun (just 
below the base of the convection zone) is of the order of 10~ 16 , 
a unit value of / would lead to Ekman layers about 10~ 8 times 
size of the domain; the numerical method used is unable to re- 
solve them. Using values of / of the order of 10 7 - 10 9 instead 
inflates the Ekman layers artificially to 10~ 4 times the size of 
the domain or larger, which can then be fully resolved. As an 
added bonus, varying / provides an easy way of varying the 
effective viscosity and thermal conductivity without changing 
the Prandtl number. As a result, the estimated values of k\ and 
k2 are unchanged from the solar value (since &2 only depends 
on the Prandtl number), while £3 and £4 are a factor of f 1 ' 2 
smaller. As long as / -C 10 15 , the hierarchy £2 <C k\ <C \k^ \,\k^\ 
is respected. 

The numerical method of solution is based on the expansion 
of the governing equations onto a spherical coordinate system 
(r, 9, 4>), followed by their projection onto Chebishev polynomi- 
als T„(co&9), and finally, solution of the resulting ODE system 
in r using a Newton-Raphson-Kantorovich algorithm. For more 
detail, see Garaud (2001). The computational domain is limited 
to the region of the radiative zone within r £ [0.02,0.7]r Q . The 
solution is found to be reasonably insensitive to the position of 
the lower boundary for this range of parameters. 

In order to study the dependence of the amplitude and depth 
of penetration of the flows on the Ekman number (i.e. on £3 
and fct), we perform a series of simulations with / ranging from 
10 7 to 10 9 . In order to verify the dependence of the amplitude 
of the solutions on £2, we perform similar calculations with a 
hypothetical Sun where the thermal conductivity in the system 
of equations (l42l is uniformly multiplied by a factor of four 
throughout the interior, resulting in a Prandtl number artificially 
decreased by a factor of four compared with its solar value. In 
that case, ki is two times smaller than in the case of the solar 



value of the Prandtl number. 

In order to study in detail the effects of the boundary con- 
ditions on the solutions, we study separately three cases where 
forcing is respectively in the azimuthal, latitudinal and radial 
directions only. In the case of radial forcing, we consider the 
two boundary conditions studied in $2] namely df / dr is null 
on the boundary, and ndf /dr = Tu r . 

The details of the expressions for the boundary conditions in 
each case are given below. 

4.2. Boundary conditions. 

In all of the simulations performed, the lower boundary con- 
ditions at r = 0.02r Q are the following: 

• impermeable condition on the radial velocity, 

• stress-free conditions on the tangential velocities ug and 

u<t>, 

• conducting condition on the temperature: we assume 
that the domain within 0.02r Q is a conducting solid 
sphere, so that the temperature perturbations satisfy 
V 2 r = within, and are regular at r = 0. We solve this 
equation with the requirement that f — > as r — > 0, and 
derive a matching condition with the temperature fluc- 
tuations at the interface. 

Near the upper boundary (at r = O.7r ), we consider the fol- 
lowing cases: 

• azimuthal forcing only: in this case, we set u r = ug = 
at the boundary (assuming no-slip conditions). Fol- 
lowing the results from helioseismology, we set u^ = 
Q.7rQsm9(l with Cl = O eq (l - acos 2 9 - Z?cos 4 9) - fi, 
where 0= fl eq (l - a/5 - 3b/35) (Gilman, Morrow & 
DeLuca, 1989), a = b = 0.15 and ft eq = 2.9 x KTV 1 . 
The temperature boundary condition is df /dr = 0. 

• latitudinal forcing only: in this case we set u r = u^ = 
at the boundary, and ug = Vosin#cos# with Vb = 1 m/s. 
Note that u^ = is guaranteed by setting a = b = 0, in 
which case the background angular velocity is ft = f2 eq = 
2.9 x 10~ 6 s _1 . The temperature boundary condition is 
df/dr = Q. 

• radial forcing only: in this case we set ug = u t j } = 0al the 
boundary and u r = Uo(l -3 cos 2 9) with Uq = lcm/s. This 
expression satisfies the global conservation of mass. 
In this last case, two temperature boundary conditions 
are explored as in $2] either df /dr = (type 1) or 
ndf /dr= Tu r (type 2). 

Note that since the governing equations are linear, the ampli- 
tudes of each component of the flow defined by (Uq, Vq,Wo) at 
the boundary can be chosen arbitrarily. Here, they were selected 
as what may be plausible flow velocities in the lower regions of 
the convection zone. 

4.3. Results 

This section summarizes the numerical results for the various 
simulations performed. 

In order to compare quantitatively the numerical solutions 
with the Cartesian model predictions we study the variation of 
the latitudinal component of the meridional flow ug, both near 
and far from the boundary, at a fixed latitude fairly close to the 
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FIG. 1 . — Latitudinal veloc ity at a latitude of 80° as a function of re-scaled depth below the convection zone £ = (r -OWtq)/ 'd-g, for the four types of boundary 
conditions described in 94.21 In each of the four plots, the solid lines (and solid symbols) correspond to simulations with a solar value of the Prandtl number, while 
the dotted lines (and open symbols) correspond to simulations with a Prandtl number artificially reduced by a factor of four. The symbols identify the value of / used 
in the simulations: a circle corresponds to / = 10 9 , a diamond to / = 10 s and a triangle to / = 10 7 . For comparison, the forcing velocities at the latitude considered 
are about u cz (80°) = 5200cm/s, v cz (80°) = 17cm/s and w cz (80°) = 0.9 cm/s. 



polar regions. We select a high latitude of about 80° since the 
Cartesian model applies best to systems in which the rotation 
axis and gravity are nearly aligned. 

In Figure [TJ we first focus on the region close to the upper 
boundary in order to single out the Ekman solution. Each plot 
corresponds to one of the four types of boundary conditions 
studied, and shows ug(r) at a latitude of about 80°. For clarity, 
the depth below the convection zone is rescaled with the Ek- 
man depth <Je'. thus, in each of the plots £ = (r— 0.7r Q )/iiE- The 
results for the solar value of the Prandtl number are shown in 
solid lines, and those corresponding to the lower value of the 
Prandtl number (equivalently, a value of k that is four times 
solar) are shown in the dotted lines. In each case, three runs 
are presented with / = 10 7 , 10 s , and 10 9 respectively and can 
be identified with the symbols. It is immediately obvious from 
observing all of these plots that there is indeed a component 
of the solution which decays and oscillates rapidly with depth 
below the convection zone on an Ekman lengthscale. From 
the conclusions of 33.61 we expect the amplitude of this rapidly 
decaying solution to scale as f/Pr in the case of azimuthal 
forcing, to be of order of the forcing velocity for any value of 
/ or Pr in the case of the latitudinal forcing, to scale as 1 / \JJ 
for radial forcing with type 1 boundary conditions and finally, 
to scale as (i/e^2 _ 1)/ Vf f° r radial forcing with type 2 bound- 
ary conditions. Note that for solar values of the background 
state and of the Prandtl number, using the standard model of 
Christensen-Dalsgaard et al. (1991), Hek 2 ~ 0.8; for a Prandtl 
number that is reduced by a factor of four, then //e^2 — 0-4. 

The comparison between these predicted scalings and the 
outcome of the numerical simulations is shown in Figure [2] 
Here, the symbols represent the amplitude of the solution de- 
fined as the maximum value achieved by |wg(r)| at the selected 
latitude, for each of the simulations performed. The lines show 
the predicted scalings, using as a reference the amplitude of the 
solution for solar values of the Prandtl number and with / = 10 7 . 
Thus, for example, the solid line in the case of azimuthal forc- 
ing only is generated by the equation 



A(f,Pr Q )=A(10 1 ,Pr G )(-L 



1/2 



(43) 



is the numerically calculated reference amplitude of the solu- 
tion for / = 10 7 and solar values of the Prandtl number. The 
dotted line in the same panel is easily calculated as 

1/2 



A(fA).25Pr.) = 2A(\(y.Pr.i(^ f 



(44) 



where A(/,Pr Q ) is the predicted amplitude of the other simula- 
tions with solar values of the Prandtl number, and A(IQ 7 ,Pr & ) 



The solid and dotted lines in the three other panels are con- 
structed in a similar fashion. It is quite clear that the predicted 
scalings fit the numerical solutions very well, except perhaps 
in the case of radial forcing with type 2 boundary conditions 
where the fit is only good to within a factor of order unity. 
The origin of this discrepancy is not entirely clear, but can be 
partly traced back to non-ideal effects in the equation of state 
(which affects the determination of H@) and to geometrical ef- 
fects (which influence the value of ki). 

In order to compare the scalings of the slowly decaying com- 
ponent of the flow with the numerical results, we now move to 
Figure[3]and Figure [4] which focus on the behavior of the solu- 
tion far from the boundary. Figure [3] shows |wg(r)| at a latitude 
of about 80° throughout the interior on a log-linear scale. The 
information about the direction of the flow (poleward or equa- 
torward) is lost in this plot, with sign reversals appearing as 
cusps in the curves pointing towards — oo. Note that it is possi- 
ble to discern the presence of the Ekman layer close to the outer 
boundary for the larger values of /. 

A quick glance at the solutions in the deep radiative interior 
reveals the behavior suggested by the Cartesian solutions: dif- 
ferent assumptions concerning the boundary conditions yield 
very different predicted flow velocities. In particular, it can 
be seen that in the case of radial forcing with type 2 boundary 
conditions (i.e. ndf /dr= Tu r ) the flows velocities are more- 
or-less independent of /, or in other words retain significant 
amplitudes for any value of the background viscosity. 

More quantitatively, from the conclusions of 33.61 we ex- 
pect the amplitude of the slowly decaying solution to scale 
as f/yPr in the case of azimuthal and latitudinal forcing, 
to scale as y/f/Pr for radial forcing with df /dr = and as 
(Hq /R)yPr in the case of radial forcing with the type 2 temper- 
ature boundary conditions. The comparison between the pre- 
dicted scalings and the outcome of the numerical simulations is 
shown in Figure [4] using the same method as described earlier 
in the case of the rapidly decaying component of the solution. 
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FIG. 2. — Comparison of the scaling of the amplitude of the rapidly decaying solution with the model predictions for the same four types of boundary conditions. 
For each simulation performed (i.e. for each type of boundary condition, for each value of / and of the Prandtl number considered), the symbols show the maximum 
value of | u$ (r)\ found numerically at the latitude considered; note that the same coding is used as in Figure[T] The solid and dotted lines shows the predicted scaling 
of this amplitude as a function of / for solar values of the Prandtl number and a quarter-times solar value of the Prandtl number respectively. The amplitude of the 
solution for solar values of the Prandtl number with / = 10 7 is used as the reference amplitude. 



Azimuthal forcing Latitudinal forcing Radial forcing, type 1 Radial forcing, type 2 




FIG. 3. — Latitudinal velocity of flows below the convection zone as a function of radius for the four types of boundary conditions described in 34.21 The line and 
symbol coding are the same as in Figure [T] 
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FIG. 4. — Comparison of the scaling of the amplitude of the slowly decaying solution with the model predictions for the same four types of boundary conditions. 
For each simulation performed (i.e. for each type of boundary condition, for each value of / and of the Prandtl number considered), the symbols show the value of 
|ng(O.5r0) found numerically at the latitude considered; again, the same coding is used. The solid and dotted lines shows the predicted scaling of this amplitude 
as a function of / for solar values of the Prandtl number and a quarter-times solar value of the Prandtl number respectively. The amplitude of the solution for solar 
values of the Prandtl number with / = 10 7 is used as the reference amplitude. 
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We can see that, as before, the predicted amplitudes agree very 
well with the numerical solutions, except perhaps in the case of 
radial forcing with type 2 boundary conditions where a small 
discrepancy remains. 

We therefore conclude that despite the simplified nature of 
the analytical analysis performed in $2] and $3] the results ob- 
tained robustly predict the scalings of the numerical solutions 
in full spherical geometry. 

Finally, in order to provide better insight into the actual solu- 
tions, we show in Figure|5]the global structure of the numerical 
results. In the left-hand-side of Figure|5J we show results for the 
solar value of the Prandtl number (and therefore solar values of 
Ic2 near the convective-radiative interface), while in the right- 
hand-side we show results for a quarter-times solar value of the 
Prandtl number. All the plots were generated for a fixed value 
of / = 10 8 (and therefore 10 8 -times solar values of the viscosity 
corresponding to 10~ 4 -times solar values of and £4). The left 
quadrants show contour-lines of the stream-function, or in other 
words, streamlines of the flow while the right quadrants show 
the angular velocity profile. The four aforementioned bound- 
ary conditions are explored: azimuthal forcing only, latitudinal 
forcing only, radial forcing only with type 1 boundary condi- 
tions (df/dr = 0) and radial forcing only with type 2 boundary 
conditions (k8T /dr= Tu r ). 

These numerical results illustrate graphically the behavior 
predicted by the Cartesian model. For example, the streamlines 
plotted all correspond to the same contours of the streamfunc- 
tion, or in other words the same mass flux. It is therefore ob- 
vious that the flows driven through latitudinal forcing are neg- 
ligible compared with the three other cases. The flows driven 
by azimuthal forcing appear to be quite strong, although this 
is merely related to the fact that the linear velocity Uq corre- 
sponding to the imposed (solar) differential rotation is more 
than two orders of magnitude larger than Vb at the same lati- 
tudes. The flows driven by radial forcing are quite strong, in 
particular given that the driving velocity Wo is at least one or- 
der of magnitude lower than Vb at the same latitude. One can 
also readily see that the flow velocities predicted using type 2 
thermal boundary conditions are much stronger than in the case 
of type 1 thermal boundary conditions. The right-side quad- 
rants reveal the effect of these meridional flows on angular mo- 
mentum transport in the interior and show the angular velocity 
profile corresponding to the ultimate steady-state achieved in 
the radiative interior should the Sun be left to evolve for many 
flow turnover times under the same applied boundary condi- 
tions. Since the turnover time can be considerably longer than 
the age of the Sun in most cases, these steady states would never 
in practise be achieved. However, they do reveal two points of 
particular interest. Firstly, they illustrate how sensitive the inte- 
rior rotation rate is to the assumed or calculated flow structure, 
and therefore also to the assumed convective-radiative interfa- 
cial boundary conditions. Secondly, they reveal the possibil- 
ity of unphysical counter-rotation permitted by the linearization 
of the problem (see 33.61 for a discussion of this effect). This 
should be taken as a cautionary warning for any linear study of 
flows and thermal perturbations in the solar radiative interior to 
always perform a self-consistency check of the validity of the 
linearization. 

5. DISCUSSION AND CONCLUSION 

We have studied the penetration of global-scale meridional 
flows generated in the solar convection zone into the radiative 



zone below. 

We have attacked the problem via a Cartesian model, as oth- 
ers have done before us, and then verified our results with more 
complete numerical modeling in the correct geometry. We have 
extended previous work done by GM04 and Mclntyre (2007) by 
taking into account the full structure of the governing equations 
and by examining the effect of boundary conditions on the solu- 
tions, or in other words, allowing a greater range of convective 
flows to act as sources for the flows in the radiative interior. 

Within a linear formalism, we confirm that the flow pattern 
in the radiative interior is a linear combination of two types of 
solutions: (a) a solution that decays away from the convective- 
radiative interface on a short length-scale related to the Ekman 
depth, and (b) a solution that decays away from the interface 
on a much longer length-scale associated with the stratification 
of the radiative interior. Our more complete treatment of the 
problem provides accurate expressions for the slowly varying 
solution. In addition, we find that the amplitudes of the rapidly 
and slowly varying solutions depend sensitively on the choice 
of dynamical and thermal interfacial forcing. 

Forcing by azimuthal and latitudinal shear at the upper 
boundary lead to both rapidly varying and slowly varying so- 
lutions, but with amplitudes that are negligible outside a few 
Ekman lengths. Forcing by direct radial pumping however 
can generate flows which are significant outside the Ekman 
boundary layer and indeed are sufficiently slowly-varying that 
they may maintain significant flow amplitudes across the whole 
depth of the radiative zone. The selected thermal boundary con- 
ditions at the convective-radiative interface are also found to 
have significant impact on the interior solution. 

In this linear problem, full solutions can be built from com- 
binations of the individual solutions obtained for each different 
type of boundary forcing that we studied. Thus in principle we 
should be able to determine the complete flow structure within 
the radiative interior. However, the problem lies in the selection 
of the interfacial conditions. We know from helioseismology 
the amplitude of the average azimuthal flows existing near the 
base of the convection zone, but know very little of the nature of 
the meridional flows or of the thermal conditions at that inter- 
face. Thus, without a complete model of the whole solar inte- 
rior which includes the turbulent solar convection zone we can 
only speculate upon the thermal and dynamical nature of the 
convective-radiative interface. Instead, when studying a broad 
range of plausible boundary conditions we find that significant 
penetration depths can be realized. 

Previous work by GM04 focused on the effects of a latitudi- 
nal source flow only. Our result is in accordance with theirs for 
this particular type of boundary condition, namely that the flows 
within the radiative interior are limited to an Ekman depth. 
However, GM04 conclude that penetration is always quenched. 
Given our analytical and numerical results in the case of other 
applied boundary conditions, we feel that this conclusion is too 
strong. Other plausible models of the interface lead to differ- 
ent results concerning the flow amplitudes within the interior. 
Without better knowledge of the conditions at the convective- 
radiative interface, it is hard to predict precisely what does hap- 
pen. 

While realistic three-dimensional global models of the so- 
lar interior combining radiative and convective regions are not 
yet numerically achievable, turbulent closure models have been 
applied to model both azimuthal and meridional flows in most 
of the solar interior (Kitchatinov & Ruediger, 2005; Rempel, 
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FIG. 5. — Global structure of the flow solutions in the radiative interior of the Sun subjected to a series of different boundary conditions at the convective-radiative 
interface (at r = 0.1 tq), for fixed value of / = 10 s . The left-hand-side set of plots shows solutions for a solar value of the Prandtl number, while the right-hand-side 
set of plots shows solutions for the case where the thermal conductivity is artificially increased by a factor of four. In each sets of plots, the streamlines are shown 
in the left-quadrant; the solid lines show clockwise flows, while the dotted lines show counter-clockwise flows. The relative angular velocity normalized by the 
equatorial value at the convective-radiative interface is shown in the right-quadrant. The "white" region in the final set of boundary conditions actually corresponds 
to an unphysical counter-rotation. This is an artefact of the linearization of the problem, combined with the high flow velocities in the interior, which would not 
occur should the fully nonlinear equations be considered. 
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2005). Using such a model, Ruediger, Kitchatinov & Arlt 
(2005) studied the penetration of meridional flows into the 
tachocline. By defining the penetration depth D pen as the dis- 
tance from the base of the convection zone to the location of 
the first reversal in ug(r), they find that D pen ~ \fEZ and con- 
clude by agreeing with GM04 that the depth of penetration of 
the flows in the interior is indeed limited to the Ekman solu- 
tion. However, it is clear from our analysis that D pen as defined 
by Ruediger, Kitchatinov & Arlt only measures the variation of 
the rapidly decaying solution (which indeed must vary on an 
Ekman length). For example, inspection of FigureQ]shows that 
at a first glance, all of the solutions appear to behave as Ek- 
man solutions, while it is only by looking more closely at the 
flows deep in the interior (in Figure [3]) that one can identify the 
presence of the slowly decaying solution with a relatively low 
but nonetheless significant amplitude. We suggest that the pre- 
dictions of the closure models for flow velocities in the interior 
could be revisited in the light of our analysis. 

Finally, the governing equations in this work have been sim- 
plified using two major approximations, which must now be 
briefly addressed. Firstly, this study is assumes the flows to be 



in a quasi-steady state. As it was shown by Mclntyre (2007) this 
assumption filters out any transient flows, which could in some 
circumstances be of much larger amplitude than these steady- 
state flows (see Spiegel & Zahn, 1992 for example). Therefore 
our conclusions are likely to underestimate the actual turnover 
time of transient meridional flows. 

Secondly, the momentum and thermal energy equation have 
been linearized with the consequences described in jj2.71 an ex- 
ample of which (the unphysical counter-rotation) is shown in 
Figure [5] We do emphasize that in the real Sun nonlinear ef- 
fects would play a role in limiting the amplitude of the merid- 
ional and azimuthal flows penetrating into the radiative zone. 
Therefore our results are consistent with the standard theory 
that no flows with turnover times faster than the thermal diffu- 
sion time are allowed into the radiative zone. This limits the 
flow turnover time in the tachocline region to a few times 10 s 
years, and into the deep interior to a few times 10 7 years. 

Neither of these caveats, however, change the main conclu- 
sion of our analysis, which is that the penetration of meridional 
flows into the radiative interior is not necessarily limited to a 
shallow Ekman depth. 
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